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Abstract 

The  Taylor-Chen-Kuszmaul  model,  which  regards  the  dynamic  fracture  process  of  brittle  materials  as  a  continuous 
accrual  of  damage,  has  been  successfully  applied  to  simulate  rock  blast  and  concrete  penetration.  This  paper  employs 
the  TCK  damage  model  to  numerical  study  on  the  effect  of  perforated  concrete  defense  layer  on  evolution  mechanism 
of  blast-induced  stress  waves.  The  numerical  results  reveal  that  the  tensile  damage  near  free  boundary  should  be  notewor¬ 
thy  under  a  higher  blast  loading,  and  the  peak  values  of  hydrostatic  pressure  beneath  an  artificial  cavity  are  largely 
reduced.  The  effects  of  cavity  dimensions  and  position  on  wave  evolution  and  reduction  are  detailedly  explored.  One 
empirical  formula  is  proposed  to  relate  the  decay  factor  of  peak  hydrostatic  pressure  to  the  dimensions  and  relative  posi¬ 
tion  of  cavity. 
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1.  Introduction 

Concrete  has  been  widely  used  for  military  and 
civil  purposes  over  many  years.  A  review  of  previ¬ 
ous  work  reveals  that  many  studies  of  concrete  shel¬ 
ters  from  the  early  1940s.  However,  most  of  this 
work  ceased  shortly  after  World  War  II  and  was 
not  resumed  until  the  1960s.  The  resistance  of  con¬ 
crete  against  dynamic  loadings  is  of  great  interest  to 
both  the  designers  of  defense  structures  and  the 
developers  of  weapon  systems  [1-3]. 

Many  findings  have  shown  that  short-duration 
high-magnitude  blast  loading  conditions  signifi- 
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cantly  influence  defense  structural  response.  It  well 
known  that  blast  loadings  are  quite  different  from 
the  conventional  loads.  For  Example,  they  are  typ¬ 
ically  applied  to  structures  at  rates  approximately 
1000  times  faster  than  earthquake-induced  loadings 
[4].  Consequently,  the  analysis  of  shelters  made  of 
brittle  concrete  subjected  to  blast  loading  becomes 
a  very  complex  issue. 

Generally  speaking,  the  blast  waves  would  evolve 
into  the  propagation  of  stress-waves  in  certain  med¬ 
ium.  The  waves  can  be  diffracted  when  hitting 
obstacles  such  as  cavity  and  crevice,  and  its  stress 
intensity  may  be  greatly  reduced  after  these  obsta¬ 
cles.  This  reduction  course  is  often  termed  as  atten¬ 
uation  of  stress-waves  [5,6]. 

New  stress-wave  migration  approaches  are  in 
expectation  because  the  attenuation  of  stress-waves 
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is  an  ongoing  topic  of  interest  in  the  design  of  con¬ 
crete  defense  structures.  A  better  understanding  of 
structural  details  behavior  under  blast  loadings  is 
very  important,  and  a  strong  interest  exists  in  devel¬ 
oping  methods  to  insure  superior  protective  perfor¬ 
mance.  For  this  sake,  some  experimental  tests 
related  have  been  performed  by  the  authors,  and  a 
series  of  field  observations  were  conducted.  Fig.  1 
is  the  appearance  of  full-scale  test,  and  Fig.  2  depicts 
the  schematic  of  the  experimental  defense  system 
conducted  in  this  project.  Fig.  3  shows  the  defense 


roof  permeated  by  blasting  cracks  after  testing. 
Detailed  presentations  of  the  full-scale  tests  were 
demonstrated  in  [7].  It  is  clear  that  a  layered  defense 
system  consists  of  three  layers:  a  soil  cover  layer,  a 
protection  layer  and  a  support  layer.  Furthermore, 
the  protection  layer  has  two  sub-layers:  a  projectile 
shelter  layer  and  a  stress  distribution  layer.  The 
stress  distribution  layer  is  inbuilt  with  artificial  cav¬ 
ity  here,  whose  function  is  to  redistribute  the  blast 
loading  over  a  larger  area.  However,  studies  on 
the  function  of  the  stress  distribution  layer  are 


Fig.  1 .  Appearance  of  field  test  for  defense  structure. 
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Fig.  2.  Schematic  of  the  defense  structure  for  one  field  test. 
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Fig.  3.  Roof  of  the  protection  layer  in  field  test. 


limited  so  far  [7-9].  Besides,  post-test  observations 
show  that  spallation  and  even  crater  are  formed  at 
the  top  of  cavities  or  near  the  non-transmission 
boundaries  of  defense  system  under  a  higher  blast 
loading.  The  formation  of  the  spallations  is  due 
mainly  to  dynamic  fracture  caused  by  the  reflected 
tensile  waves.  Those  phenomena  ignored  in  the  pre¬ 
vious  studies  [7,9]  will  be  paid  more  attention  to  in 
this  work. 

The  remainder  of  this  paper  is  organized  as  fol¬ 
lows:  In  Section  2,  the  TCK  continuum  damage 
model  [10,11]  for  brittle  concrete  is  introduced  and 
succinctly  implemented  into  LS-DYNA  [12].  The 
prerequisites  for  numerical  analyses  are  given  in  Sec¬ 
tion  3,  including  the  verification  of  the  newly-imple¬ 
mented  numerical  tool.  Section  4  numerically 
explores  the  effect  of  a  single  rectangular  cavity  on 
the  attenuation  of  stress-waves,  and  one  empirical 
formula  is  brought  forward  to  relate  the  decay  factor 
for  peak  P  to  the  cavity  parameters.  At  last,  some 
conclusions  and  remarks  are  stated  in  Section  5. 

The  primary  purpose  of  this  paper  is  to  investi¬ 
gate  the  dynamic  responses  of  perforated  concrete 
defense  systems  using  the  TCK  damage  model  and 
provide  references  to  the  designers  of  intelligent 
protective  shelters. 

2.  Description  of  TCK  damage  model 

Concrete  is  typically  brittle  and  non-homoge- 
neous  in  composition.  The  basic  assumption  of  the 
damage  model  is  that  the  concrete  is  permeated  by 
an  array  of  randomly  distributed  cracks  which  grow 


and  interact  with  one  another  under  tensile  loading. 
The  damage  is  reflected  in  the  degradation  of  the 
material  stiffness  following  the  equations  derived 
in  [13]  for  a  random  array  of  penny-shaped  cracks 
in  an  isotropic  elastic  medium.  The  effective  bulk 
modulus  of  the  cracked  solid  is 

17  i  c 

-=l--/(»)Cd  (l) 

where  f(v)  =  (1  —  f2)/(l  —  2v),  K  and  v  are  the  ori¬ 
ginal  elastic  constants  (bulk  modulus  and  Possion’s 
ratio)  for  intact  material,  and  barred  quantities 
represent  the  corresponding  degraded  constants 
for  the  damaged  material.  Cd  is  the  crack  density 
parameter,  which  is  assumed  to  be  proportional  to 
the  product  of  N ,  the  number  of  cracks  of  per  unit 
volume,  and  r3,  the  cube  of  the  average  crack 
dimension  in  a  representative  volume.  Following 
the  work  in  [14],  N  is  expressed  as  a  Weibull  statis¬ 
tical  distribution  function  activated  by  the  bulk 
strain  measure  sv  =  (sx  +  sy  +  sz)/3,  according  to 

N  =  k(sv)m  (2) 

where  k  and  m  are  material  constants  to  be  deter¬ 
mined  from  strain-rate  dependent  tensile  fracture 
stress  data. 

To  extent  the  range  of  crack  densities,  an 
expression 

V  =  vexp  0<  ,8  <  1  (3) 
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based  on  the  percolation  theory  [15]  was  used  for 
the  equivalent  Poisson’s  ratio  v .  The  value  of  /?  con¬ 
trols  the  unloading  and  reloading  behaviors  and  re¬ 
laxes  the  restriction  of  elastic  unloading  in  the 
original  model  [16]. 

Based  on  kinetic  energy  consideration,  an  expres¬ 
sion  is  derived  [17]  for  the  nominal  fragment  radius 
for  dynamic  fragmentation  in  a  brittle  material  as: 

(4) 

2  \^pC£v max  J 

In  Eq.  (4)  p  is  the  mass  density,  C  is  the  uniaxial 
wave  speed  ^E/p^j,  svmSLX  is  the  maximum  volu¬ 
metric  strain  rate  experienced  by  the  representative 
volume  element  for  the  entire  loading  process,  and 
Kic  is  the  fracture  toughness  of  the  material. 

When  bulk  tension  occurs  in  the  material,  the 
crack  density  parameter  Cd  can  be  calculated  as: 

Cd  =  2^)m(jf)  e„L  (5) 

From  Eq.  (1),  the  damage  scalar  can  be  calcu¬ 
lated  as  follows: 

D  =  jfAv)C'd  (6) 

The  degraded  shear  modulus  G  and  Bulk  modu¬ 
lus  K  are  written  as: 


G=G{  1  -D) 
K  =  K(  1  -D) 


(7) 


Using  Eq.  (7),  the  elastic  strain  relations  for 
stress  are  then  adopted  in  K  —  G  form  for  bulk  ten¬ 
sion  of  concrete: 


A  Gij  =  KAskkdij  +  2GAeij  (8) 

where  Sy  =  0  for  /  ^  j  and  Sy  =  1  for  i  =  j.  Aoy  and 
Ae-y  are  increment  for  the  stress  tensor  and  the  devi- 
atoric  strain  tensor,  respectively,  and  Avkk  is  the  vol¬ 
umetric  strain  increment. 

The  plastic/kinematic  algorithm  used  here  for  the 
response  of  concrete  in  compression  is  the  usual 
elastic  radial-return  approach  found  in  most  nonlin¬ 
ear  finite  element  codes  [18]. 


The  von  Mises  yield  condition  is  expressed  as: 


where  a  «  y  |  (sy  —  &>*/)(%  —  (%),  oy  denotes  the 
yield  stress  which  is  assumed  constant. 


The  center  of  the  yield  in  deviatoric  stress  space  is 
given  by  tensor  w  which  has  initial  value  of  zero. 
The  increment  of  o  is  given  by: 

Aa>ij  =  ^£pef/k  (10) 

where  d t  is  the  time  step,  Ep  is  the  plastic  hardening 
modulus,  and  the  plastic  strain  ratio  £?•  is  the  differ¬ 
ence  between  total  and  elastic  strain  (right  super¬ 
script  e)  ratios: 


Table  1  gives  the  main  concrete  parameters  for 
this  concrete  damage  model  [11]. 

3.  Preconditions  for  numerical  simulations 

3.1.  Equivalent  blast  loading 

As  indicated  in  Fig.  4,  the  planar  charge  in  the 
field  research  program  is  adopted.  Because  the  rise 
time  of  blast  loading  is  very  short  and  the  overpres¬ 
sure  decays  exponentially  after  the  peak  value,  an 
equivalent  triangular  pressure  pulse  (see  Fig.  5)  is 
used  to  approximate  the  blast  loading  [9,19,20].  In 
the  present  study,  the  peak  value  (9max  is  taken  as 
120.0  MPa  at  time  tmax  =  0.0  (is  (micro  second) 
and  then  decreases  linearly  to  zero  at  time 
tc  =  100.0  (is.  It  is  in  response  to  the  display  of  pres¬ 
sure  sensors  S5  and  S6  (see  Fig.  2)  and  directly 
applied  on  the  upper  surface  of  the  projectile  shelter 
layer. 

3.2.  Implementation  of  the  constitutive  model 

The  hydrocode  FS-DYNA  from  the  Livermore 
Software  Technology  Corporation  is  used  for  the 
numerical  analyses  in  this  article.  It  can  handle 
either  two-  or  three-  dimensional  problems  for 
explosion  and  high-speed  impact  [12,21]. 


Table  1 

Material  parameters  for  concrete 


P  (kg.m  3) 

E  (GPa) 

V 

<ry  (MPa) 

k  (m"3) 

M 

^IC  (MPa.m-1/2) 

2520.0 

20.68 

0.18 

20.0 

5.753e21 

6.0 

2.747 
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Fig.  4.  Field  assembly  of  blasting  fuses  for  planar  charge. 


The  TCK  model  described  above  has  been  imple¬ 
mented  as  a  constitutive  augmentation  to  the 
dynamic  finite  code  LS-DYNA,  and  the  detailed 
flowchart  is  illustrated  in  Fig.  6.  The  implementa¬ 
tion  follows  general  guidelines  for  such  augmenta¬ 
tion  and  employs  the  given  interface  subroutines. 
The  augmentation  is  vectorized  consistent  with  the 
basic  LS-DYNA  code,  and  incorporates  damage 
as  a  history  variable. 


shows  a  comparison  between  the  numerical  predic¬ 
tions  with  the  measured  data  of  radial  stress  in  place 
at  a  radial  distance  of  3.0  m  and  a  depth  of  2.5  m.  It 
can  be  found  that  the  developed  user  subroutine  is 
capable  of  calculating  the  stress  with  good  accuracy. 

Fig.  2  shows  the  schematic  drawing  for  one 
field  experiment.  A  series  of  on-site  observations 
were  conducted.  In  numerical  simulations,  only  the 
right-half  domain  is  considered,  and  a  total  of  3636 
quadrilateral  elements  are  used,  including  180  ele¬ 
ments  for  steel  plate.  The  steel  plate  is  modeled  with 
material  Type  3  available  in  LS-DYNA.  It  is  a  bi-lin- 
ear  elastoplastic  model.  Table  2  lists  the  material 
parameters  for  the  steel  plate.  The  blast  loading 
shown  in  Fig.  5  is  applied  upon  the  projectile  shelter 
layer.  A  comparison  between  the  measured  with  cal¬ 
culated  wave  forms  of  hydrostatic  pressure  in  the 
place  of  sensor  S8  (see  Fig.  2)  is  shown  in  Fig.  8. 
The  sensor  S8  is  at  a  depth  of  1.25  m  beneath  the 
lower  side  of  cavity.  It  can  be  seen  that  the  numerical 
simulation  of  the  peak  value  and  waveform  show  a 
consistent  comparison  with  the  measured  data. 
However,  the  tensile  phenomenon  is  obvious  in  the 
curve  tail  for  the  numerical  case.  The  reason  is  that 
the  lower  boundary  of  the  calculation  domain  is 
regarded  as  a  reflected  one  (i.e.  non-transmission) 
in  simulations.  In  fact,  there  exist  soil  layers  beneath 
the  defense  system. 

4.  Attenuation  effect  of  single  cavity  on  stress-waves 

In  the  field  tests,  multiple  rectangular  cavities 
were  adopted  (refer  to  Fig.  2).  How  about  the  effect 
of  single  cavity  on  attenuation  of  stress-waves?  This 
section  will  investigate  it  based  on  the  above  verified 
numerical  tool. 

4.1.  Definition  of  decay  factor 


3.3.  Verification  of  user  define  subroutine 

This  model  was  ever  used  in  a  study  of  blast 
induced  damage  in  rock  by  Taylor  et  al.  [10].  The 
single  column  of  high  explosive  (Iregel  1175U)  was 
bottom  detonated.  It  was  2.5  m  in  height,  0.081  m 
in  radius,  and  was  centered  at  3.75  m  below  the 
ground  surface.  The  oil  shale  was  assumed  to  be  iso¬ 
tropic  and  homogenous  material  with  pre-existing 
micro  cracks.  Included  in  [10]  were  the  constitutive 
parameters  for  the  rock  and  the  explosive.  Fig.  7 


To  evaluate  the  attenuation  extent  of  a  stress- 
wave  due  to  a  cavity  in  a  concrete  defense  system, 
a  dimensionless  variable,  ‘decay  factor’  DF,  is 
defined  as  follows: 


where  o0  denotes  the  peak  value  of  stress  compo¬ 
nent,  such  as  t tx ,  t xy,  and  hydrostatic  pressure  P  at 
a  specified  position  without  cavity  in  the  defense 
layer,  and  o  is  the  peak  value  of  the  same  compo¬ 
nent  at  the  same  position  with  a  cavity  in  the 
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Fig.  6.  Implementation  flowchart  for  the  TCK  model. 


Fig.  7.  Comparison  between  calculated  and  measured  data  for 
radial  stress. 


Fig.  8.  Comparison  between  calculated  and  measured  data  for 
hydrostatic  pressure. 


Table  2 


Material  parameters  for  steel  plate 

P  (kg  m  3) 

E  (GPa) 

V 

ay  (MPa) 

7900.0 

207.0 

0.29 

340.0 

defense  layer.  Eq.  (12)  reveals  that  a  larger  DF 
would  imply  a  better  attenuation  effect. 


4.2.  Computational  domain  with  cavity 

Fig.  9  shows  a  typical  concrete  defense  layer 
embedded  with  a  rectangular  cavity.  The  computa¬ 
tional  domain  for  this  plane  strain  model  is  taken 
as  600  cm  x  300  cm  (width  x  height).  The  embedded 
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Fig.  9.  Simulation  of  a  perforated  stress  distribution  layer. 


cavity  cross-section  has  the  dimensions  of  21  in 
length  and  h  in  width.  Other  dimensions  are  indi¬ 
cated  in  Fig.  9.  Introducing  the  term  relative  dis¬ 
tance  y  as  the  distance  between  point  A  and  the 
lower  side  of  the  cavity,  the  dimension  Fhas  the  fol¬ 
lowing  relationship: 

Y  =y  +  e  +  h  (13) 

A  half  domain  is  considered  here  owing  to  sym¬ 
metry.  As  an  example,  the  right-half  domain  shown 
in  Fig.  9  has  the  following  boundary  conditions: 
traction  boundary  on  the  top,  non-transmission 
boundaries  on  the  right  and  lower  sides,  and  sym¬ 
metric  boundary  on  the  left  side. 

4.3.  Analysis  of  cavity  parameters  on  stress-waves 

4.3.1.  Effect  of  cavity  length  l 

Five  lengths  are  taken  in  the  computations:  0.00, 
0.30,  0.60,  0.90  and  1.20  m.  Other  geometrical 
parameters  are  fixed  as  e  =  0.6  m  and  h  —  0.3  m. 
Fig.  10  shows  the  effect  of  different  lengths  /  on 
the  attenuation  of  the  peak  P.  For  comparisons, 
the  computations  are  also  carried  out  based  on  the 
linear-elastic  (L-E)  model.  The  variations  of  peak 
P  with  relative  distance  y  (refer  to  definition  shown 
in  Fig.  9)  are  shown  in  Fig.  10(a)  for  the  L-E  model 
and  Fig.  10(b)  for  the  TCK  model.  The  results  show 
that  the  peak  P  increases  with  relative  distance  y 
and  decreases  with  increasing  length  /.  It  can  also 
be  seen  that  the  elastoplastic  properties  of  concrete 
plays  a  crucial  role  on  the  magnitude  of  peak  P. 


Relative  distance  y  (m) 

(a)  L-E  model 


Relative  distance  y  (m) 

(b)  TCK  model 


Fig.  10.  Effect  of  different  lengths  /  of  cavity  on  peak  hydrostatic 
pressure. 
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For  instance,  when  /  =  0.30  m  and  y  =  0.80  m,  the 
peak  P  is  24.25  MPa  for  the  L-E  model  and  only 
9.29  MPa  for  the  TCK  model.  The  corresponding 
relationship  between  DF  for  the  peak  P  and  relative 
distance  y  is  shown  in  Fig.  1 1 .  It  can  be  established 
that  DF  decreases  with  increasing  y  and  increases 
with  cavity  length  for  the  same  y. 

Fig.  12  shows  the  hydrostatic  pressure  contours 
for  the  cases  with  or  without  cavity.  A  higher  ten¬ 
sile  stress  region  is  observed  above  the  cavity,  and 
the  maximum  tension  is  about  12.1  MPa  at 
t  =  347.8  (is.  The  peak  P  beneath  the  cavity  is  evi¬ 
dently  lower  and  even  approaches  zero.  In  addition, 
the  reflected  tensile  waves  can  also  be  found  near  the 
left  and  right  boundaries.  As  for  the  non-cavity  case 
shown  in  Fig.  12(b),  it  is  evident  that  the  continuity 
of  wave  front  is  hardly  interrupted  and  the  tensile 
stress  zones  only  rise  at  the  side  boundaries  of  the 
domain. 

Assuming  the  duration  of  blast  loading  propor¬ 
tional  to  its  peak  value,  Fig.  13  compares  the  simu¬ 
lated  tensile  damage  for  different  situations.  Clearly, 


Relative  distance  y  (m) 


Fig.  11.  Variation  of  decay  factors  for  peak  P  with  different 
lengths  of  cavity  for  TCK  model. 


more  serious  tensile  damage  above  the  cavity  and 
near  the  boundaries  can  be  noticed  with  the  incre¬ 
ment  of  blast  loading.  It  is  because  the  reflected  ten- 
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Fig.  12.  Contours  of  hydrostatic  pressure  for  cases  with  and  without  cavity  at  the  same  time  (Unit:  Mbar.  IMbar  =  105  MPa). 


Z.L.  Wang,  Y.C.  Li  /  Theoretical  and  Applied  Fracture  Mechanics  47  (2007)  15-25 


23 


TCK  DAMAGE  MODEL 
Time  =  506.05 


Fringe  Levels 
1.000e+00 
9.000e  01 
8.000e  01 
7.000e  01 
6.000e-01 
5.000e-01 
4.000e-01 
3.000e  01 
2.000e-01 
I.OOOe  01 
0.000e+00 


1 


(a)  pmM=120MPa,  rc=100Ms 


TCK  DAMAGE  MODEL 
Time  =  504.09 


(b)  Qmm  =240MPa,  tc=  200  ms 


Fringe  Levels 
1.000e+00  _ 
9.000e-01  _ 
8.000e  01  _ 
7.000e-01  _ 
G.OOOe-OI  _ 
5.000e-01  _ 
4.000e-01  _ 
3.000e-01  _ 
2.000e-01  _ 
I.OOOe  01  J 
O.OOOetOO  I 


TCK  DAMAGE  MODEL 
Time  =  502 


(c)  0max=36OMPa,  /c=300Ms 


Fringe  Levels 
1.000e+00_ 
9.000e-01  _ 
8.000e  01  _ 
7.000e  01  _ 
6.000e-01  _ 
5.000e  01  _ 
4.000e  01  _ 
3.000e-01  _ 
2.000e-01  _ 
l.OOOe-OI  _ 
0.000e+00 J 


Fig.  13.  Development  of  tensile  damage  with  blast  loading  (/  =  0.60  m,  h  =  0.30  m  and  e  =  0.60  m)  at  t  =  500  ps. 


sile  waves  become  stronger  under  the  higher  exter¬ 
nal  loadings. 

4.3.2.  Effect  of  cavity  width  h 

Four  values  of  h  =  0.30,  0.60,  0.90,  1.20  m,  when 
l=e  =  0.30  m,  are  considered  to  study  the  effect  of 
cavity  width  on  the  attenuation  of  stress- waves.  The 
numerical  results  are  shown  in  Fig.  14.  It  is  noted 


that  the  cavity  width  h  has  little  influence  on  the 
decay  factor  for  the  peak  P. 

4.3.3.  Effect  of  cavity  position  e 

The  effect  of  cavity  position  e  is  examined  here. 
The  e  values  are  taken  as  0.30,  0.60,  0.90  and 
1.20  m  when  l  =  h  =  0.30  m.  Fig.  15  compares  the 
curves  of  P  versus  t.  For  all  e  values,  the  peak  values 
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Relative  distance  y  (m) 


Fig.  14.  Variation  of  decay  factors  for  peak  P  with  different 
widths  h  of  cavity. 


Time  (j.is) 


Fig.  15.  Effect  of  cavity  position  e  on  hydrostatic  pressure  history 
at  y  =  0.5  m  point  on  F-axis. 

do  not  vary  much,  but  the  arrival  time  of  the  peak 
values  gradually  lags  with  an  increase  of  the  propa¬ 
gation  distance  of  stress-waves.  Fig.  16  shows  the 
decay  factor  of  the  peak  oy  versus  relative  distance 
y.  It  is  evident  that  all  the  curves  are  almost  identi¬ 
cal.  Therefore,  the  position  e  of  the  cavity  has  little 
influence  on  the  decay  factor  for  the  peak  oy. 

4.3.4.  Empirical  formula  for  decay  factor 

Based  on  the  above  numerical  results,  the  effects 
of  h  and  e  on  the  peak  P  are  established  to  be  insig¬ 
nificant.  Thus,  DF  can  be  considered  as  a  function 
of  length  /  and  relative  distance  y: 

DF  =  f(y,  l)  (14) 


Relative  distance  y  (m) 


Fig.  16.  Variation  of  decay  factors  for  peak  P  with  different 
positions  e  of  cavity. 

In  order  to  determine  the  particular  form  of  Eq. 
(14),  one  value  of  e  (0.30  m),  four  values  of  /  (0.30, 
0.60,  0.90  and  1.20  m)  and  five  values  of  y  (0.2,  0.4, 
0.6,  0.8  and  1.0  m)  are  chosen  from  Fig.  11. 
Through  a  curve  fitting  regression  analysis,  the 
DF  for  the  peak  P  at  a  specified  point  on  the  F-axis 
is  approximately  expressed  as  follows: 

DF  =  exp(-|7)  (15) 

This  formula  shows  that  DF  increases  with  /  and 
decreases  with  y.  When  term  y  approaches  zero,  DF 
is  almost  equal  to  unity.  Fig.  17  shows  a  comparison 
between  the  numerical  results  and  the  predictions 
from  the  empirical  formula.  Good  agreement  can 
be  observed.  Therefore,  Eq.  (15)  can  be  used  to  esti¬ 
mate  the  length  /  of  a  rectangular  cavity  if  the 


Fig.  17.  Comparison  of  empirical  formulae  and  numerical  results 
of  decay  factors  for  peak  P. 
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expected  DF  at  a  relative  distance  y  is  given.  When 
DF  and  /  are  given,  Eq.  (15)  can  also  be  used  to  esti¬ 
mate  the  nearest  safe  distance  (i.e.  jvn). 

5.  Conclusions 

Continuous  demands  on  safety  and  intelligent 
design  have  placed  increasing  emphasis  on  perfor¬ 
mance  analyses  of  concrete  defense  structures.  This 
paper  consists  of  three  aspects:  user-defined  consti¬ 
tutive  subroutine,  numerical  investigations  and 
empirical  formula.  From  these  studies,  the  following 
conclusions  and  understandings  may  be  drawn. 

(1)  The  TCK  model  is  realized  by  coding  it  as  user 
subroutine  in  the  dynamic  finite  element  pro¬ 
gram  LS-DYNA.  It  not  only  can  capture  the 
tensile  damage  above  cavities  or  near  non¬ 
transmission  boundaries  but  also  can  evaluate 
the  elastoplastic  behaviors  of  concrete  in 
compression. 

(2)  Artificial  cavities  embedded  in  the  concrete 
defense  layer  exert  significant  influences  on 
the  attenuation  of  stress-waves.  For  a  rectan¬ 
gular  cavity,  the  decay  factor  largely  depends 
on  its  length,  while  its  width  has  little  effect 
if  the  width  is  larger  than  0.30  m  in  our 
computations. 

(3)  One  empirical  formula  is  obtained  to  correlate 
the  decay  factor  for  peak  P,  dimensions  of  a 
cavity  and  relative  distance.  The  formula  can 
be  employed  to  determine  the  minimum  length 
of  a  cavity  /m in,  or  minimum  safe  distance  jVin 
when  the  decay  factor  is  given.  Therefore,  it 
can  serve  as  useful  tools  in  the  design  of  intel¬ 
ligent  defense  engineering. 

Acknowledgements 

This  work  reported  herein  was  partially  sup¬ 
ported  by  the  Postdoctoral  Science  Foundation  of 
China  (Project  no.2004036468).  The  authors  also 
wish  to  express  their  gratitude  to  The  Third  Institute 
of  the  Corps  of  Engineers  for  providing  materials  of 
filed  experiments. 

References 

[1]  R.P.  Kenedey,  A  review  of  procedure  for  the  analysis  and 
design  of  concrete  structure  to  resisting  missile  impact  effects, 
Nuclear  Engineering  37  (1976)  183-203. 


[2]  D.V.  Balandin,  N.N.  Bolotnik,  W.D.  Pilkey,  Optimal 
Protection  from  Impact,  Shock,  and  Vibration,  Gordon 
and  Breach  Science  Publisher,  Toronto,  2001. 

[3]  T.L.  Teng,  Y.A.  Chu,  F.A.  Chang,  H.S.  Chin,  Numerical 
analysis  of  oblique  impact  on  reinforced  concrete,  Cement  & 
Concrete  Composites  27  (2005)  481-492. 

[4]  T.  Krauthammer,  Blast-resistant  structural  concrete  and 
steel  connections,  International  Journal  of  Impact  Engineer¬ 
ing  22  (1999)  887-910. 

[5]  T.P.  Philippidis,  T.G.  Aggelis,  Experimental  study  of  wave 
dispersion  and  attenuation  in  concrete,  Ultrasonics  43  (2005) 
584-595. 

[6]  Z.L.  Wang,  Y.C.  Li,  J.G.  Wang,  Study  of  stress  waves  in 
geomedia  and  effect  of  a  soil  cover  layer  on  wave  attenuation 
using  a  ID  finite  difference  method,  Computers  and  Geo¬ 
sciences  32  (2006)  1535-1543. 

[7]  Y.C.  Li,  X.J.  Wang,  X.Z.  Hu,  Study  on  layered  design  and 
its  defending  function,  Technical  Report  of  National 
Defense  Engineering,  USTC,  China,  2004. 

[8]  B.  Rohani,  Shielding  methodology  for  conventional  kinetic 
energy  weapons,  Technical  Report  SL-8f-8.  U.S.  Army 
engineers  Corps,  Waterways  Experimental  Station,  Vicks¬ 
burg.  M.S,  1987. 

[9]  Z.L.  Wang,  J.G.  Wang,  Y.C.  Li,  C.F.  Leung,  Attenuation 
effect  of  artificial  cavity  on  air-blast  waves  in  an  intelligent 
defense  layer,  Computers  and  Geotechnics  33  (2006)  132-141. 

[10]  L.M.  Taylor,  E.P.  Chen,  J.S.  Kuszmaul,  Micro-crack 
induced  damage  accumulation  in  brittle  rock  under  dynamic 
loading,  Computer  Methods  in  Applied  Mechanics  and 
Engineering  55  (1986)  301-320. 

[11]  E.P.  Chen,  Simulation  of  concrete  perforation  based  on  a 
continuum  damage  model,  in:  A.  Carpinteri  (Ed.),  Size-Scale 
Effects  in  the  Failure  Mechanisms  of  Materials  and  Struc¬ 
tures,  E  &  FN  Spon  Press,  1996. 

[12]  Livemore  Software  Technology  Corporation  (LSTC),  LS- 
DYNA  Keyword  User’s  manual,  California,  USA,  2003. 

[13]  B.  Budiansky,  R.J.  O’Connell,  Elastic  moduli  of  a  cracked 
solid,  International  Journal  of  Solids  and  Structures  12 
(1976)  81-97. 

[14]  D.E.  Grady,  M.E.  Kipp,  Continuum  modelling  of  explosive 
fracture  in  oil  shale,  International  Journal  of  Rock  Mechan¬ 
ics  and  Mining  Sciences  17  (1980)  147-157. 

[15]  R.  Englman,  Z.  Jaeger,  Theoretical  aids  for  the  improvement 
of  blasting  efficiencies  in  oil  shale  and  rocks,  AP-TR-12/87, 
Soreq  Nuclear  Research  Center,  Yavne,  Israe,1987. 

[16]  E.P.  Chen,  Non-local  effects  on  dynamic  damage  accumu¬ 
lation  in  brittle  solids,  International  Journal  for  Numerical 
and  Analytical  Methods  in  Geomechanics  23  (1999)  1-21. 

[17]  D.E.  Grady,  The  mechanics  of  fracture  under  high-rate  stress 
loading:  In  Bazant,  Z.P.  (Ed.),  Preprints  of  the  William 
Prager  Symp.  on  Mechanics  of  Geomaterials:  Rocks,  Con¬ 
crete  and  Soils,  Northwestern  University,  Evanston,  IL, 
1983. 

[18]  D.R.J.  Owen,  E.  Hinton,  Finite  elements  in  plasticity: 
Theory  and  Practice,  Pineidge  Press,  Swansea,  UK,  1980. 

[19]  S.  Guruprasad,  A.  Mukherjee,  Layered  sacrificial  claddings 
under  blast  loading,  Part  I  -  analytical  studies,  International 
Journal  of  Impact  Engineering  24  (2000)  957-973. 

[20]  P.S.  Bulson,  Explosive  loading  of  engineering  structures,  E  & 
FN  Spon  Press,  London,  1997. 

[21]  J.O.  Hallquist,  LS-DYNA  Theoretical  Manual,  California, 
USA,  1998. 


